#include <cuda_runtime.h>
#include <cublas_v2.h>
#include <cusolverDn.h>
#include <stdio.h>
#include"library1.h"
#include<string.h>
#include<stdlib.h>
#include <time.h>
#include<math.h>      //gcc�����ʱ����Ҫ��-lmѡ��

#define THREASHOLD 1e-8
#define ITERATION 11
#define COL 1
#define NUM_THD_SORT 256
#include "sort.hxx"

__device__ int sign(double number) {
    if(number<0)
        return -1;
    else
        return 1;
}

int procSize;   //���������Ŀ

/*�Ѿ���д���ļ�*/
void writeToFile(double *matrix,int rows,int columns,char* file){
    FILE *fp;
    if((fp=fopen(file,"w"))==NULL){
        perror("fopen");
        exit(1);
    }
    fprintf(fp,"%d\t%d\n",rows,columns);        //�ļ����м�¼���������������
    int i,j;
    for(i=0;i<rows;++i){
        for(j=0;j<columns;++j){
            fprintf(fp,"%-10f\t",matrix[i*columns+j]);
        }
        fprintf(fp,"\n");
    }
    fclose(fp);
}

/*������д���ļ�*/
void vectorToFile(double *vector,int len,char* file){
    FILE *fp;
    if((fp=fopen(file,"w"))==NULL){
        perror("fopen");
        exit(1);
    }
    int i;
    for(i=0;i<len;++i){
        fprintf(fp,"%-10f\t",vector[i]);
    }
    fclose(fp);
}

__device__ void sort(double *a,int *p,int max){
    __shared__ double tnorm[700];
    __shared__ int tp[700];
    int tid=threadIdx.x;
    /*while(tid<max){
        printf("a[%d]=%lf,p[%d]=%d\n",tid,a[tid],tid,p[tid]);
        tid=tid+blockDim.x;
    }
    tid=threadIdx.x;*/
    int i,j,k;
    int n=1;  //�鲢��������Ԫ�ظ���
    while(n<max){
        while(tid<max){
            if((tid%(n*2)==0)){
                for(i=tid,j=tid+n,k=tid;;k++){
                    if(k>=max)
                        break;
                    //printf("a[%d]=%lf\n",tid,a[tid]);
                    if(a[i]>a[j]||j>=max){
                        /*if(p[i]>500)
                            printf("1:tid=%d,a[%d]=%lf,p[%d]=%d,tp[%d]=%d,%d\n",tid,i,a[i],i,p[i],k,tp[k],n);*/
                        tnorm[k]=a[i];
                        tp[k]=p[i];
                        i++;
                        if(i==tid+n){
                            for(;j<tid+n+n;j++){
                                k++;
                                if(k>=max)
                                    break;
                                /*if(p[i]>500)
                                    printf("2:tid=%d,a[%d]=%lf,p[%d]=%d,tp[%d]=%d,%d\n",tid,j,a[j],j,p[j],k,tp[k],n);*/
                                tnorm[k]=a[j];
                                tp[k]=p[j];
                                
                            }
                            break;
                        }
                    }else{
                        /*if(p[i]>500)
                            printf("23:tid=,%d,a[%d]=%lf,p[%d]=%d,tp[%d]=%d,%d\n",tid,j,a[j],j,p[j],k,tp[k],n);*/
                        tnorm[k]=a[j];
                        tp[k]=p[j];
                        j++;
                        if(j==tid+n+n){
                            for(;i<tid+n;i++){
                                k++;
                                if(k>=max)
                                    break;
                                /*if(p[i]>500)
                                    printf("15:tid=%d,a[%d]=%lf,p[%d]=%d,tp[%d]=%d,%d\n",tid,i,a[i],i,p[i],k,tp[k],n);*/
                                tnorm[k]=a[i];
                                tp[k]=p[i];
                            }
                            break;
                        }
                    }
                }
            }
            tid=tid+blockDim.x;
        }
        n*=2;
        __syncthreads();
        tid=threadIdx.x;
        while(tid<max){
            /*if(tp[tid]>500){
                printf("a[%d]=%lf,p[%d]=%d,tp[%d]=%d,%d\n",tid,a[tid],tid,p[tid],tid,tp[tid],threadIdx.x);
            }*/
            a[tid]=tnorm[tid];
            p[tid]=tp[tid];
            /*if(p[tid]>500){
                printf("a[%d]=%lf,p[%d]=%d,tp[%d]=%d,%d,%d\n",tid,a[tid],tid,p[tid],tid,tp[tid],n,threadIdx.x);
            }*/
            tid=tid+blockDim.x;
        }
        tid=threadIdx.x;
        __syncthreads();
    }
}

__device__ void push_forward(int *p_a,int *p_b,int i,int length_b){
    //swap last i & spin p_b
    volatile int temp;
    __shared__ int p[500];
    int tid=threadIdx.x;
    //swap last i element
    if(threadIdx.x==0){
        temp=p_a[length_b-i-1];
        p_a[length_b-i-1]=p_b[length_b-i-1];
        p_b[length_b-i-1]=temp;
        temp = p_b[0];
    }
    __syncthreads();
    //spin p_b to p
    while(tid<length_b-1){
        p[tid]=p_b[tid+1];
        tid+=blockDim.x;
    }
    __syncthreads();
    if(threadIdx.x==0)
        p[length_b-1]=temp;
    tid=threadIdx.x;
    __syncthreads();
    //save p_b
    while(tid<length_b){
        p_b[tid]=p[tid];
        tid+=blockDim.x;
    }
    __syncthreads();
}

__device__ void push_forward_odd(int *p_a,int *p_b,int length_b){
    //p_a & p_b make a whole ring
    volatile int temp;
    __shared__ int p[500];
    int tid=threadIdx.x;
    //spin p_a to p
    if(threadIdx.x==0){
        temp=p_a[length_b];
    }
    __syncthreads();
    while(tid<length_b){
        p[tid+1]=p_a[tid];
        tid+=blockDim.x;
    }
    tid=threadIdx.x;
    __syncthreads();
    //swap element 0
    if(threadIdx.x==0)
        p[0]=p_b[0];
    __syncthreads();
    //save p_a
    while(tid<=length_b){
        p_a[tid]=p[tid];
        tid+=blockDim.x;
    }
    tid=threadIdx.x;
    __syncthreads();
    //spin p_b
    while(tid<length_b-1){
        p[tid]=p_b[tid+1];
        tid+=blockDim.x;
    }
    tid=threadIdx.x;
    __syncthreads();
    if(threadIdx.x==0)
        p[length_b-1]=temp;
    __syncthreads();
    while(tid<length_b){
        p_b[tid]=p[tid];
        tid+=blockDim.x;
    }
    __syncthreads();
    }

__device__ int p_a[2][700];
__device__ int p1[700];
__device__ int p_ab[2][700][700];

//�����������ڻ�
__device__ void vectorProduct1(double *a,double *b,int len, double* cache){
    //ÿһ�����϶���cache�����Ŀ������໥֮�䲻Ӱ��
    cache[threadIdx.x]=0.0;
    //tidΪ�̵߳�ƫ����
    int tid=threadIdx.x;
    double temp=0.0;
    while(tid<len){
        temp+=a[tid]*b[tid];
        tid+=blockDim.x;
    }
    
    cache[threadIdx.x] = temp;
    __syncthreads();

    //�ڿ��ڼ��㲿�ֺ�
    if(threadIdx.x<128)
        cache[threadIdx.x]+=cache[threadIdx.x+128];
    __syncthreads();
    if(threadIdx.x<64)
        cache[threadIdx.x]+=cache[threadIdx.x+64];
    __syncthreads();
    if(threadIdx.x<32)
        cache[threadIdx.x]+=cache[threadIdx.x+32];
    if(threadIdx.x<16)
        cache[threadIdx.x]+=cache[threadIdx.x+16];
    if(threadIdx.x<8)
        cache[threadIdx.x]+=cache[threadIdx.x+8];
    if(threadIdx.x<4)
        cache[threadIdx.x]+=cache[threadIdx.x+4];
    if(threadIdx.x<2)
        cache[threadIdx.x]+=cache[threadIdx.x+2];
    if(threadIdx.x == 0);
        cache[0] = cache[0] + cache[1];
}


__global__ void sort_and_save(double *A,double *V,double *U,double *S,double *A1,double *V1,int row,int col,size_t pitch){
    __shared__ double d[700];
    __shared__ int p[700];
    int tid=threadIdx.x;
    while(tid<row){
        double* r=(double*)((char*)A+tid*pitch);
        d[tid]=sqrt(vectorProduct(r,r,row));
        p[tid]=tid;
        tid=tid+blockDim.x;
    }
    __syncthreads();
    sort(d,p,row);
    __syncthreads();
    tid=threadIdx.x;
    while(tid<row){
        S[tid]=d[tid];
        //printf("S[%d]=%lf,d[%d]=%lf,p[%d]=%d\n",threadIdx.x,S[threadIdx.x],threadIdx.x,d[threadIdx.x],threadIdx.x,p[threadIdx.x]);
        double* Ai=(double*)((char*)A1+tid*pitch);
        double* Aj=(double*)((char*)A+p[tid]*pitch);
        double* Vi=(double*)((char*)V1+tid*pitch);
        double* Vj=(double*)((char*)V+p[tid]*pitch);
        for(int i=0;i<row;i++){
            Ai[i]=Aj[i];
            Vi[i]=Vj[i];
        }
        tid=tid+blockDim.x;
    }
}

__global__ void getU(double *A1,double *S,double *U,int row,int col,int pitch){
    int tid=threadIdx.x;
    while(tid<row){
        if(S[tid]<THREASHOLD)
            for(int i=0;i<col;++i)
                U[tid*col+i]=0.0;
        else
            for(int i=0;i<col;++i){
                double* Ai=(double*)((char*)A1+tid*pitch);
                U[tid*col+i]=Ai[i]/S[tid];
            }
        tid=tid+blockDim.x;
    }
}

__device__ volatile int Arrayin[500];
__device__ void __gpu_syncNew(volatile int *Arrayin, int iter)
{
    // only thread0 is used for synchronization
    if ( threadIdx.x == 0 )
    {

        Arrayin[blockIdx.x] += 1;
        __threadfence();

        if ( blockIdx.x < gridDim.x-1 )
        {
            while ( Arrayin[blockIdx.x+1] < iter )
            {
                // Do nothing here
            }
        }
            
        if ( blockIdx.x == gridDim.x-1 )
        {
            while ( Arrayin[0] < iter )
            {
                // Do nothing here
            }
        }

    }

    __syncthreads();
}

__device__ double d1[700];
__device__ double d2[700];
__global__ void getD(double *A,int row,int col,size_t pitch){
    __shared__ double cache[512];
    double* r=(double*)((char*)A+blockIdx.x*pitch);
    vectorProduct1(r,r,row,cache);
    if(threadIdx.x==0){
        d1[blockIdx.x]=cache[0];
        d2[blockIdx.x]=cache[0];
    }
}
__global__ void getrank(int row,int col,int *pass){
    __shared__ int p[700];
    int tid=threadIdx.x;
    while(tid<row){
        pass[tid]=1;
        p[tid]=tid;
        tid=tid+blockDim.x;
    }
    __syncthreads();
    sort(d2,p,row);
    __syncthreads();

    int i = threadIdx.x;
    while(i < row){ //整个if block效率偏低，考虑重新设计
        //printf("p[%d]=%d,d[%d]=%lf\n",i,p[i],i,d[i]);
        if(i%2==0){
            p_a[0][i/2]=p[i];
            //printf("p_a[0][%d]=%d, p_ab[0][%d][0]=%d\n", i, p_a[0][a],i, p_ab[0][a][0]);
        }
        else{
            p_a[1][i/2]=p[i];
            //printf("p_a[1][%d]=%d, p_ab[1][%d][1]=%d\n", i, p_a[1][b],i, p_ab[1][b][0]);
        }

        p1[i] = p[i];
        i += blockDim.x;
    }
    __syncthreads();
}

__global__ void create_index_pair(int row, int col){
    __shared__ int p[700];
    int tid;
    //先给p赋值再用这个值是啥意思？？
    tid=threadIdx.x;
    while(tid<row){
        p[tid]=tid;
        tid=tid+blockDim.x;
    }
    __syncthreads();

    if( threadIdx.x == 0 ){
        int h=0, k=0;
        for ( int i = 0; i < row; i ++ ){
            if ( i%2 == 0 ){
                p_ab[0][h][0]=p[i];
                p_a[0][h++]=p[i];
            }
            else{
                p_ab[1][k][0]=p[i];
                p_a[1][k++]=p[i];
            }
        }
    }
    __syncthreads();

    int count = 0;
    for ( int i = 1; i < row-(1-row%2); i ++ ) {

        if (row % 2 == 0) {

            if ((i - 1 == 2 * (count + 1)) && threadIdx.x == 0) {
                count = count + 1;
            }

            __syncthreads();

            push_forward(p_a[0], p_a[1], count, row / 2);
        } else {
            push_forward_odd(p_a[0], p_a[1], row / 2);
        }

        __syncthreads();
        tid = threadIdx.x;
        while (tid < row / 2) {
            p_ab[0][tid][i] = p_a[0][tid];
            p_ab[1][tid][i] = p_a[1][tid];
            tid = tid + blockDim.x;
        }
        __syncthreads();

    }
}

__global__ void gpusvd(double *A,double *V, int row, int col,size_t pitch,int *pass,int *count_i) {       //�Ǽ��ɲ���svd�ֽ�
    __shared__ double ele, ele1, ele2;
    __shared__ double cache[256];
    __shared__ double Acol[2][500];
    __shared__ double Vcol[2][500];
    __shared__ int hold_vector;
    int pa[2];
    //int pp1[500];
    //int pp2[500];
    int tid;
    // reset flags, signal waiting blocks
    if (threadIdx.x == 0){
        Arrayin[blockIdx.x] = 0;
    }
    __syncthreads();

    double *Aco[2],*Vco[2];
    for(int i=0;i<row-(1-row%2);i++){
        //pa[0]=p_a[0][blockIdx.x];
        //pa[1]=p_a[1][blockIdx.x];
        
        // determind which pair to process
        pa[0]=p1[p_ab[0][blockIdx.x][i]];
        pa[1]=p1[p_ab[1][blockIdx.x][i]];

        Aco[0]=(double*)((char*)A+pa[0]*pitch);
        Aco[1]=(double*)((char*)A+pa[1]*pitch);
        Vco[0]=(double*)((char*)V+pa[0]*pitch);
        Vco[1]=(double*)((char*)V+pa[1]*pitch);

        // load vector into shared
        if( i==0 ){
            tid=threadIdx.x;
            while(tid<row){
                Acol[0][tid]=Aco[0][tid];
                Acol[1][tid]=Aco[1][tid];
                Vcol[0][tid]=Vco[0][tid];
                Vcol[1][tid]=Vco[1][tid];
                tid=tid+blockDim.x;
            }
            __syncthreads();
        }
        else{
            tid=threadIdx.x;
            while(tid<row){
                if(hold_vector==0){ //hold current vector to pair with another
                    Acol[0][tid]=Acol[1][tid];
                }
                Acol[1][tid]=Aco[1][tid];
                Vcol[1][tid]=Vco[1][tid];
                tid=tid+blockDim.x;
            }
            __syncthreads();
        }

        if(threadIdx.x==0){
            hold_vector=(blockIdx.x!=(row/2-1-i/2));
        }
        __syncthreads();

        vectorProduct1(Acol[0],Acol[1],row,cache);                                       //5
        if(threadIdx.x==0){
            ele=cache[0];
        }
        __syncthreads();

        if(fabs(ele)>THREASHOLD){
            if(threadIdx.x==0){
                ele1=d1[pa[0]];
                ele2=d1[pa[1]];
                pass[blockIdx.x]=0;
            }
            __syncthreads();

            double aro;
            double aco;
            double tao=(ele1-ele2)/(2*ele);
            double tan=sign(tao)/(fabs(tao)+sqrt(1+pow(tao,2)));
            double cos=1/sqrt(1+pow(tan,2));
            double sin=cos*tan;

            if ( hold_vector == 0 ){
                tid = threadIdx.x;
                while (tid < row) {
                    aro = Acol[0][tid] * cos + Acol[1][tid] * sin;
                    aco = Acol[1][tid] * cos - Acol[0][tid] * sin;
                    Acol[0][tid] = aro;
                    Acol[1][tid] = aco;
                    aro = Vcol[0][tid] * cos + Vcol[1][tid] * sin;
                    aco = Vcol[1][tid] * cos - Vcol[0][tid] * sin;
                    Vcol[0][tid] = aro;
                    Vcol[1][tid] = aco;

                    Aco[0][tid] = Acol[0][tid];
                    Vco[0][tid] = Vcol[0][tid];
                    tid = tid + blockDim.x;
                }
                __syncthreads();
            }
            else{
                tid = threadIdx.x;
                while (tid < row) {
                    aro = Acol[0][tid] * cos + Acol[1][tid] * sin;
                    aco = Acol[1][tid] * cos - Acol[0][tid] * sin;
                    Acol[0][tid] = aro;
                    Acol[1][tid] = aco;
                    aro = Vcol[0][tid] * cos + Vcol[1][tid] * sin;
                    aco = Vcol[1][tid] * cos - Vcol[0][tid] * sin;
                    Vcol[0][tid] = aro;
                    Vcol[1][tid] = aco;

                    Aco[1][tid] = Acol[1][tid];
                    Vco[1][tid] = Vcol[1][tid];
                    tid = tid + blockDim.x;
                }
                __syncthreads();
            }

            if(threadIdx.x==0){
                d1[pa[0]]=ele1+ele*tan;
                d1[pa[1]]=ele2-ele*tan;
            }
            __threadfence();
        }
        else if ( hold_vector == 0 ){
            tid = threadIdx.x;
            while (tid < row) {
                Aco[0][tid] = Acol[0][tid];
                Vco[0][tid] = Vcol[0][tid];
                tid = tid + blockDim.x;
            }
            __syncthreads();
        }

        __gpu_syncNew(Arrayin, i+1);
        //__gpu_sync(Arrayin);
    }

    if ( hold_vector == 0 ){
        tid = threadIdx.x;
        while (tid < row) {
            Aco[1][tid] = Acol[1][tid];
            Vco[1][tid] = Vcol[1][tid];
            tid = tid + blockDim.x;
        }
        __syncthreads();
    }
    else{
        tid = threadIdx.x;
        while (tid < row) {
            Aco[0][tid] = Acol[0][tid];
            Vco[0][tid] = Vcol[0][tid];
            tid = tid + blockDim.x;
        }
        __syncthreads();
    }
}

int main(int argc, char*argv[]){
    //int *dev_procSize;
    int width,height;
    /*׼�������A��V*/
    FILE *fp;
    fp=fopen(argv[1],"r");//���ļ�
    fscanf(fp,"%d",&height);
    fscanf(fp,"%d",&width);
    double *host_A=(double*)malloc(width*height*sizeof(double));
    double *host_V=(double*)malloc(height*height*sizeof(double));
    for(int i=0;i<height;i++) //������
        for(int j=0;j<height;j++)
            fscanf(fp,"%lf",&host_A[i*width+j]);
    fclose(fp);//�ر��ļ�
    clock_t t1 = clock(); //��ʼ��ʱ��
    procSize = height/(COL*2);
    double *dev_A,*dev_A1,*dev_S,*dev_U,*dev_V,*dev_V1;
    for(int j=0;j<height;++j){
        host_V[j*height+j]=1.0;
    }
    double *host_U=(double*)malloc(width*width*sizeof(double));
    double *host_S=(double*)malloc(height*sizeof(double));
    size_t pitch;
    int *pass,*dev_pass,*count,*dev_count;
    cudaEvent_t start1;  
    cudaEventCreate(&start1);  
    cudaEvent_t stop1;  
    cudaEventCreate(&stop1);
    cudaHostAlloc((void**)&pass,height*sizeof(int),cudaHostAllocWriteCombined|cudaHostAllocMapped);
    cudaHostAlloc((void**)&count,sizeof(int),cudaHostAllocMapped);
    cudaHostGetDevicePointer(&dev_pass,pass,0);
    cudaHostGetDevicePointer(&dev_count,count,0);
    cudaMallocPitch(&dev_A,&pitch,width*sizeof(double),height);
    cudaMallocPitch(&dev_V,&pitch,height*sizeof(double),height);
    cudaMallocPitch(&dev_A1,&pitch,width*sizeof(double),height);
    cudaMallocPitch(&dev_V1,&pitch,height*sizeof(double),height);
    cudaMalloc((void**)&dev_U,width*width*sizeof(double));
    cudaMalloc((void**)&dev_S,height*sizeof(double));
    cudaMemcpy2D(dev_A,pitch,host_A,width*sizeof(double),width*sizeof(double),height,cudaMemcpyHostToDevice);
    cudaMemcpy2D(dev_V,pitch,host_V,height*sizeof(double),height*sizeof(double),height,cudaMemcpyHostToDevice);
    cudaEventRecord(start1, NULL);


    create_index_pair<<<1,width/2>>>(height,width);
    getD<<<height,256>>>(dev_A,height,width,pitch);
    cudaDeviceSynchronize();
    int iteration= ITERATION;
    while(iteration-- > 0)
    {
        //printf("iter = %d\n", iteration);
        count[0]=0;
        getrank<<<1,256>>>(height,width,dev_pass);
        cudaDeviceSynchronize();
        gpusvd<<<procSize,256>>>(dev_A,dev_V,height,width,pitch,dev_pass,dev_count);
        cudaDeviceSynchronize();
        int allpass=0;
        for(int i=0;i<procSize;i++){
            allpass=allpass+pass[i];
        }

        if(allpass==procSize){
            break;
        }
    }
    sort_and_save<<<1,256>>>(dev_A,dev_V,dev_U,dev_S,dev_A1,dev_V1,height,width,pitch);
    cudaDeviceSynchronize();
    getU<<<1,256>>>(dev_A1,dev_S,dev_U,height,width,pitch);
    cudaDeviceSynchronize();
    printf("number of iteration=%d\n",ITERATION-iteration-1);
    cudaEventRecord(stop1, NULL);  
    cudaEventSynchronize(stop1);  
    float msecTotal1 = 0.0f;  
    cudaEventElapsedTime(&msecTotal1, start1, stop1); 
    printf("time = %f\n",msecTotal1);
    cudaError_t  error_check101 ;
    error_check101 = cudaGetLastError();
    if( error_check101 != cudaSuccess ){
        printf("%s\n" , cudaGetErrorString( error_check101 ) );
    }
    cudaMemcpy2D(host_V,height*sizeof(double),dev_V1,pitch,height*sizeof(double),height,cudaMemcpyDeviceToHost);
    cudaMemcpy(host_S,dev_S,height*sizeof(double),cudaMemcpyDeviceToHost);
    cudaMemcpy(host_U,dev_U,width*width*sizeof(double),cudaMemcpyDeviceToHost);
    clock_t t2 = clock(); //������ʱ
    char outFileV[7]={'V','X','.','m','a','t','\0'};
    writeToFile(host_V,height,height,outFileV);
    char outFileU[7]={'U','X','.','m','a','t','\0'};
    writeToFile(host_U,width,width,outFileU);
    char outFileS[7]={'S','X','.','m','a','t','\0'};
    vectorToFile(host_S,height,outFileS);
    printf("cputime = %lf\n",(double)((t2-t1)*1000/CLOCKS_PER_SEC));
    cudaFree(dev_A);
    cudaFree(dev_V);
    cudaFree(dev_A1);
    cudaFree(dev_V1);
    cudaFree(dev_U);
    cudaFree(dev_S);
    cudaFree((void*)count);
    cudaFree((void*)pass);
    free(host_A);
    free(host_V);
    free(host_U);
    free(host_S);
    return 0;
}
